---
title: "7-3 All Cause Mortality Day 28"
description: |
Analyse of the all-cause mortality secondary outcome in ASCOT.
author:
- name: Rob Mahar
affiliation: University of Melbourne
- name: James Totterdell
affiliation: University of Sydney
date: today
freeze: true
---
# Preamble {#preamble}
```{r}
#| label: pkgs
#| code-summary: Load packages
library (tidyverse)
library (labelled)
library (kableExtra)
library (cmdstanr)
library (posterior)
library (bayestestR)
library (bayesplot)
library (matrixStats)
library (lubridate)
library (patchwork)
library (ggdist)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: data
#| code-summary: Prepare dataset
devtools:: load_all ()
all_dat <- read_all_no_daily ()
# Anticoagulation set (intention-to-treat, including missing)
acs_itt_dat <- all_dat %>%
filter_acs_itt () %>%
transmute_model_cols_grp_aus_nz ()
# Anticoagulation set (intention-to-treat, excluding missing)
acs_itt_nona_dat <- acs_itt_dat %>%
filter (! is.na (D28_death))
# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
filter_concurrent_intermediate ()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>%
filter (! is.na (D28_death))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
filter_concurrent_std_aspirin ()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>%
filter (! is.na (D28_death))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
filter_concurrent_therapeutic ()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>%
filter (! is.na (D28_death))
```
```{r}
#| label: models
#| code-summary: Load the Stan models
model_full <- cmdstan_model (file.path ("stan" , "binary" , "logistic_site_epoch.stan" ))
model_site_only <- cmdstan_model (file.path ("stan" , "binary" , "logistic_site.stan" ))
model_epoch_only <- cmdstan_model (file.path ("stan" , "binary" , "logistic_epoch.stan" ))
model_fixed_only <- cmdstan_model (file.path ("stan" , "binary" , "logistic.stan" ))
```
```{r}
#| label: helpers
#| code-summary: Helper functions for this analysis script
make_stan_data <- function (
dat,
vars = NULL ,
outcome = NULL ,
beta_sd = NULL ,
ctr = contr.orthonorm) {
X <- make_X_design (dat, vars, ctr)
nXassign <- sum (grepl ("rand" , colnames (X))) - 1
beta_sd <- c (2.5 , rep (1 , nXassign), beta_sd)
y <- dat[[outcome]]
N <- dim (X)[1 ]
K <- dim (X)[2 ]
epoch <- dat$ epoch
M_epoch <- max (dat$ epoch)
region <- dat[["ctry_num" ]]
M_region <- max (region)
site <- dat[["site_num" ]]
M_site <- max (site)
region_by_site <- region_by_site <- dat %>%
dplyr:: count (ctry_num, site_num) %>%
pull (ctry_num)
list (N = N, K = K, X = X, y = y,
M_region = M_region, region = region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
beta_sd = beta_sd)
}
make_stan_data_concurrent <- function (
dat,
vars = NULL ,
outcome = NULL ,
beta_sd = NULL ,
ctr = contr.orthonorm) {
X <- make_X_design (dat, vars, ctr, includeA = FALSE )
nXassign <- sum (grepl ("rand" , colnames (X))) - 1
beta_sd <- c (2.5 , rep (1 , nXassign), beta_sd)
y <- dat[[outcome]]
N <- dim (X)[1 ]
K <- dim (X)[2 ]
list (N = N, K = K, X = X, y = y, beta_sd = beta_sd)
}
make_mortality_table <- function (dat, format = "html" ) {
tab <- dat %>%
mutate (CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C3" , "C4" ),
labels = str_replace (intervention_labels ()$ CAssignment[- 1 ], "<br>" , " " ))) %>%
group_by (CAssignment) %>%
summarise (
Patients = n (),
Known = sprintf (
"%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
Missing = sprintf (
"%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Died by day 28 ` = sprintf (
"%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sprintf (
"%i (%.1f)" , sum (! is.na (D28_death)), 100 * mean (! is.na (D28_death))),
Missing = sprintf (
"%i (%.1f)" , sum (is.na (D28_death)), 100 * mean (is.na (D28_death))),
` Died by day 28 ` = sprintf (
"%i (%.1f)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
)
) %>%
mutate (CAssignment = fct_inorder (CAssignment)) %>%
gather (key, value, - CAssignment, factor_key = T) %>%
spread (key, value)
colnames (tab)[1 ] <- "n ( \\ %)"
if (format == "latex" ) {
colnames (tab) <- linebreak (colnames (tab), align = "c" , linebreaker = "<br>" )
}
kable (tab,
format = format,
align = "lrrrrr" ,
escape = FALSE ,
booktabs = TRUE
) %>%
kable_styling (font_size = 10 , latex_options = "HOLD_position" ) %>%
row_spec (nrow (tab), bold = TRUE )
}
```
# Descriptive
## Anticoagulation
```{r}
#| code-summary: Summary table
#| label: descriptive-table-mortality-anticoagulation
#| tbl-cap: Summary of mortality to day 28 by anticoagulation intervention.
save_tex_table (make_mortality_table (
all_dat %>%
filter_acs_itt (),
"latex" ),
"outcomes/secondary/anticoagulation-summary-mortality" )
make_mortality_table (
all_dat %>%
filter_acs_itt ()
)
```
## Age
```{r}
#| label: fig-age-po
#| code-summary: Relationship age to outcome
#| fig-cap: |
#| Relationship (logistic regression linear in age)
#| between age at entry and the primary outcome.
agedat <- acs_itt_dat %>%
dplyr:: count (D28_death, AgeAtEntry) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 0 ` + ` 1 ` + ` <NA> ` ,
p = ` 1 ` / (` 1 ` + ` 0 ` ))
agemod <- glm (
cbind (` 1 ` , ` 0 ` ) ~ AgeAtEntry,
data = agedat,
family = binomial ())
agedat <- agedat %>%
mutate (
ypred = predict (agemod, newdata = agedat, type = "response" )
)
p1 <- ggplot (agedat, aes (AgeAtEntry, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" )
p2 <- ggplot (agedat, aes (AgeAtEntry, p)) +
geom_point () +
geom_vline (xintercept = 60 , linetype = 2 ) +
geom_line (aes (y = ypred)) +
labs (y = "Proportion \n died by day 28" , x = "Age at entry" )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-age.pdf" ), p, height = 2 , width = 6 )
p
```
## Country
```{r}
#| label: fig-country-7-3
#| code-summary: Relationship country to outcome
#| fig-cap: Mortality by country.
tdat <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (Country = fct_infreq (PT_CountryName), D28_death) %>%
group_by (Country) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
)
p1 <- ggplot (tdat, aes (Country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (tdat, aes (Country, p_1)) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Country of enrolment" ) +
ylim (0 , NA )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-country.pdf" ), p, height = 2 , width = 6 )
p
```
## Site
```{r}
#| label: fig-site-7-3
#| code-summary: Relationship site to outcome
#| fig-cap: Day 28 mortality by site within country.
tdat <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (
Country_lab = Country,
Site_lab = fct_infreq (Location),
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" )),
Site = PT_LocationName,
D28_death) %>%
group_by (Country, Site) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (
n = ` 1 ` + ` 0 ` + ` <NA> ` ,
p_1 = ` 1 ` / (` 1 ` + ` 0 ` ),
p_miss = ` <NA> ` / (` 1 ` + ` 0 ` + ` <NA> ` )
) %>%
ungroup ()
p1 <- ggplot (tdat, aes (Site_lab, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (tdat, aes (Site_lab, p_1)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Site of enrolment" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p <- p1 / p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-cal-po
#| code-summary: Relationship calendar date to outcome
#| fig-cap: |
#| Relationship between calendar date and the primary outcome.
caldat <- all_dat %>%
filter_acs_itt () %>%
dplyr:: count (D28_death, yr = year (RandDate), mth = month (RandDate)) %>%
spread (D28_death, n, fill = 0 ) %>%
mutate (p = ` 1 ` / (` 1 ` + ` 0 ` ),
n = ` 1 ` + ` 0 ` + ` <NA> ` )
p1 <- ggplot (caldat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_point () +
labs (
y = "Proportion \n died by day 28" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 ) +
ylim (0 , 0.15 )
p2 <- ggplot (caldat, aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p <- p2 | p1
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-3-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
# Modelling
The analysis for all cause mortality is specified equivalently to that for the primary outcome. It includes intervention as randomised, age category, country, site nested within country, epoch, and intervention eligibility. The main analysis restricts the population to only those participants randomised to the anticoagulation domain.
## Primary model
```{r}
#| label: model1-sampling
#| code-summary: Data and sampling
seed <- 59579
mdat <- make_stan_data (
dat = acs_itt_nona_dat,
vars = c ("inelgc3" , "agegte60" , "ctry" ),
outcome = "D28_death" ,
beta_sd = c (10 , 2.5 , 1 , 1 ))
snk <- capture.output (
mfit <- model_full$ sample (
data = mdat,
seed = seed,
refresh = 0 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
chains = 8 ,
parallel_chains = 8 ,
adapt_delta = 0.99
))
mfit$ save_object (file.path ("outputs" , "models" , "secondary" , "7-3.rds" ))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (mdat$ X)
names (mdrws$ gamma_epoch) <- acs_itt_nona_dat %>% dplyr:: count (epoch, epoch_lab) %>% pull (epoch_lab)
names (mdrws$ gamma_site) <- acs_itt_nona_dat %>% dplyr:: count (site_num, site) %>% pull (site)
site_facet <- factor (mdat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
mdrws$ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]
mdrws$ Atrt <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c (
"Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),
"Therapeutic vs low" = exp (mdrws$ Ctrt[3 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[3 ]),
"Low with aspirin vs therapeutic" = exp (mdrws$ Ctrt[2 ] - mdrws$ Ctrt[3 ])
)
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
setNames (exp (mdrws$ beta[7 : 8 ]), c ("Ineligible aspirin" , "Age \u2265 60" )),
setNames (exp (mdrws$ beta[9 : 10 ]), c ("Australia/New Zealand" , "Nepal" ))
)
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save
#| code-summary: Odds ratio summary table
save_tex_table (
odds_ratio_summary_table (mdrws$ OR, "latex" ),
"outcomes/secondary/7-3-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
#| label: odds-ratio-summary-primary-model-epoch-site-terms
#| code-summary: Epoch and site terms
p <- plot_epoch_site_terms (
exp (mdrws$ gamma_epoch),
exp (mdrws$ gamma_site),
site_facet
)
pth <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" , "7-3-primary-model-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
```{r}
#| fig-cap: Posterior densities for treatment comparisons.
p <- plot_or_densities (mdrws$ compare)
p
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" , "gamma_site" , "gamma_epoch" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["gamma_site" ])
mcmc_trace (mdrws["gamma_epoch" ])
mcmc_trace (mdrws[c ("tau_site" , "tau_epoch" )])
```
:::
### Posterior predictive
```{r}
#| label: primary-acs-itt-ppc
#| fig-cap: Posterior predictive distribution versus observed proportion (red diamond).
y_ppc <- mdrws$ y_ppc
ppc_dat <- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc))
grp_ppc <- function (grp) {
ppc_dat %>%
group_by (grp = !! grp) %>%
summarise (
mean_y = mean (D28_death),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
}
plot_grp_ppc <- function (dat, lab = "" ) {
dat %>%
ggplot (aes (y = grp, xdist = rvar_mean_y_ppc)) +
stat_interval (size = 2 ) +
geom_point (aes (x = mean_y, y = 1 : nrow (dat)), colour = "red" , shape = 23 ) +
labs (
x = "Posterior predictive \n proportion" ,
y = lab)
}
ppc_C <- grp_ppc (sym ("CAssignment" ))
ppc_ctry <- grp_ppc (sym ("country" ))
ppc_epoch <- grp_ppc (sym ("epoch" ))
ppc_site <- ppc_dat %>%
group_by (Country = country, grp = site_raw) %>%
summarise (
mean_y = mean (D28_death),
rvar_mean_y_ppc = rvar_mean (y_ppc)
)
p1 <- plot_grp_ppc (ppc_C, "Anticoagulation \n intervention" ) + labs (x = "" )
p2 <- plot_grp_ppc (ppc_ctry, "Country" ) + labs (x = "" )
p3 <- plot_grp_ppc (ppc_epoch, "Epoch" ) + labs (x = "" )
p4 <- plot_grp_ppc (ppc_site %>% filter (Country == "IN" ), "Sites \n India" ) + labs (x = "" )
p5 <- plot_grp_ppc (ppc_site %>% filter (Country == "AU" ), "Sites \n Australia" ) + labs (x = "" )
p6 <- plot_grp_ppc (ppc_site %>% filter (Country == "NP" ), "Sites \n Nepal" )
p7 <- plot_grp_ppc (ppc_site %>% filter (Country == "NZ" ), "Sites \n NZ" )
p <- ((p3 | p1 / p2) /
( (p4 | p5) / (p6 | p7) +
plot_layout (heights = c (3 , 1 )) ) ) +
plot_layout (heights = c (1 , 1.5 ), guides = "collect" )
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" ,
"7-3-primary-model-acs-itt-ppc.pdf" )
ggsave (pth, p, width = 6 , height = 5.5 )
p
```
## Sensitivity: Reduced model (no site terms)
The primary model is re-fit without the hierarchical site terms.
```{r}
#| code-summary: Fit model
mdat <- make_stan_data (
acs_itt_nona_dat,
c ("inelgc3" , "agegte60" , "ctry" ),
"D28_death" ,
c (10 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- model_epoch_only$ sample (
data = mdat,
seed = 14646 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.99
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" , "gamma_epoch" , "tau_epoch" )))
names (mdrws$ beta) <- colnames (mdat$ X)
epoch_map <- acs_itt_nona_dat %>% dplyr:: count (epoch, epoch_lab)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 4 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
setNames (exp (mdrws$ beta[- (1 : 6 )]), c ("Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for no epoch model.
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
plot_epoch_terms (exp (mdrws$ gamma_epoch))
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" , "gamma_epoch" , "tau_epoch" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["gamma_epoch" ])
mcmc_trace (mdrws["tau_epoch" ])
```
:::
## Sensitivity: Reduced model (no epoch terms)
The primary model is re-fit without the epoch terms.
```{r}
#| code-summary: Fit model
mdat <- make_stan_data (
acs_itt_nona_dat,
c ("inelgc3" , "agegte60" , "ctry" ),
"D28_death" ,
c (10 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- model_site_only$ sample (
data = mdat,
seed = 123513 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.99
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" , "gamma_site" , "tau_site" )))
site_map <- acs_itt_nona_dat %>% dplyr:: count (site_num, site)
site_facet <- factor (mdat$ region_by_site, labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
names (mdrws$ beta) <- colnames (mdat$ X)
names (mdrws$ gamma_site) <- site_map$ site
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 4 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
setNames (exp (mdrws$ beta[- (1 : 6 )]), c ("Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for no epoch model.
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
plot_site_terms (exp (mdrws$ gamma_site), site_facet)
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" , "gamma_site" , "tau_site" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["gamma_site" ])
mcmc_trace (mdrws["tau_site" ])
```
:::
## Sensitivity: Reduced model (no site or epoch)
The primary model is re-fit without the site and epoch terms.
```{r}
#| code-summary: Fit reduced model (no site or epoch)
mdat <- make_stan_data (
acs_itt_nona_dat,
c ("inelgc3" , "agegte60" , "ctry" ),
"D28_death" ,
c (10 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- model_fixed_only$ sample (
data = mdat,
seed = 12415 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.99
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" )))
names (mdrws$ beta) <- colnames (mdat$ X)
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 4 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
setNames (exp (mdrws$ beta[- (1 : 6 )]), c ("Ineligible aspirin" , "Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for no epoch model.
odds_ratio_summary_table (mdrws$ OR)
```
:::{.callout-caution collapse="true"}
### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
```
:::
## Sensitivity: Concurrent Controls
A sensitivity analysis is performed on the reduced concurrent control analysis sets for each intervention.
### Intermediate
- **Set**: ACS-ITT-intermediate
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Overview of outcome
#| tbl-cap: Summary of outcome for intermediate dose concurrent enrolments.
save_tex_table (
make_mortality_table (acs_itt_concurc2_dat, "latex" ),
"outcomes/secondary/7-3-anticoagulation-concurrent-intermediate-summary" )
make_mortality_table (acs_itt_concurc2_dat)
```
```{r}
#| code-summary: Fit model
mdat <- make_stan_data_concurrent (
acs_itt_concurc2_nona_dat,
c ("agegte60" , "ctry" ),
"D28_death" ,
c (2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- model_fixed_only$ sample (
data = mdat,
seed = 38135 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.95
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" )))
names (mdrws$ beta) <- colnames (mdat$ X)
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]))
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" )),
setNames (exp (mdrws$ beta[- (1 : 2 )]), c ("Age \u2265 60" , "Australia/New Zealand" , "Nepal" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for model.
save_tex_table (
odds_ratio_summary_table (mdrws$ OR, "latex" ),
"outcomes/secondary/7-3-primary-model-acs-itt-concurrent-intermediate-summary-table" )
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-intermediate
p <- plot_or_densities (mdrws$ compare)
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
```
:::
### Low-dose with aspirin
- **Set**: ACS-ITT-aspirin
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Overview of outcome
#| tbl-cap: Summary of outcome for low-dose with aspirin concurrent enrolments.
save_tex_table (
make_mortality_table (acs_itt_concurc3_dat, "latex" ),
"outcomes/secondary/7-3-anticoagulation-concurrent-stdaspirin-summary" )
make_mortality_table (acs_itt_concurc3_dat)
```
```{r}
#| code-summary: Fit model
mdat <- make_stan_data_concurrent (
acs_itt_concurc3_nona_dat,
c ("agegte60" , "ctry" ),
"D28_death" ,
c (2.5 , 1 )
)
snk <- capture.output (
mfit <- model_fixed_only$ sample (
data = mdat,
seed = 13578 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" )))
names (mdrws$ beta) <- colnames (mdat$ X)
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 3 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Low with aspirin" )),
setNames (exp (mdrws$ beta[- (1 : 3 )]), c ("Age \u2265 60" , "Australia/New Zealand" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for model.
save_tex_table (
odds_ratio_summary_table (mdrws$ OR, "latex" ),
"outcomes/secondary/7-3-primary-model-acs-itt-concurrent-stdaspirin-summary-table" )
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-aspirin
p <- plot_or_densities (mdrws$ compare)
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
```
:::
### Therapeutic
- **Set**: ACS-ITT-therapeutic
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Overview of outcome
#| tbl-cap: Summary of outcome for therapeutic dose concurrent enrolments.
save_tex_table (
make_mortality_table (acs_itt_concurc4_dat, "latex" ),
"outcomes/secondary/7-3-anticoagulation-concurrent-therapeutic-summary" )
make_mortality_table (acs_itt_concurc4_dat)
```
```{r}
#| code-summary: Fit model
mdat <- make_stan_data_concurrent (
acs_itt_concurc4_nona_dat,
c ("agegte60" , "ctry" ),
"D28_death" ,
c (2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- model_fixed_only$ sample (
data = mdat,
seed = 13578 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws (c ("beta" )))
names (mdrws$ beta) <- colnames (mdat$ X)
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[2 : 3 ]
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Therapeutic vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c (
setNames (exp (mdrws$ Ctrt), c ("Intermediate" , "Therapeutic" )),
setNames (exp (mdrws$ beta[- (1 : 3 )]), c ("Age \u2265 60" , "India" , "Australia/New Zealand" ))
)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Odds ratio summary table for model.
save_tex_table (
odds_ratio_summary_table (mdrws$ OR, "latex" ),
"outcomes/secondary/7-3-primary-model-acs-itt-concurrent-therapeutic-summary-table" )
odds_ratio_summary_table (mdrws$ OR)
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Posterior densities for treatment comparisons, ACS-ITT-therapeutic
p <- plot_or_densities (mdrws$ compare)
p
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (variables = c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
```
:::